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1. Introduction 

The modular polynomial, 4'e{x, y), parameterizes pairs of elliptic curves with a cyclic isogeny of 
degree i between them. Modular polynomials provide the defining equations for modular curves, 
and are useful in many different aspects of computational number theory and cryptography. For 
example, computations with modular polynomials have been used to speed elliptic curve point- 
counting algorithms f |BSS99j Chapter VII). 

The standard method for computing modular polynomials consists of computing the Fourier ex- 
pansion of the modular j-function and solving a linear system of equations to obtain the integral 
coefficients of (^^{x^y). According to Elkies (see |Elk98j §3) this method has a running time of 
0(£^+'^) if one uses fast multiplication. However, our analysis (given in the Appendix) shows that 
the running time of this method is, in fact, 6(^9/2+^) using fast multiplication. 

The object of the current paper is to compute the modular polynomial (for prime t) directly mod- 
ulo a prime p, without first computing the coefficients as integers. Once the modular polynomial 
has been computed for enough small primes, our approach can also be combined with the Chinese 
Remainder Theorem (CRT) approach as in |CJNST98| or |ALVn3j to obtain the modular polyno- 
mial with integral coefficients or with coefficients modulo a much larger prime using Explicit CRT. 
Our algorithm does not involve computing Fourier coefficients of modular functions. The running 
time of our algorithm turns out to be 0{l^^'') using fast multiplication. We believe our method 
is interesting as it is asymptotically faster; and is an essentially different approach to computing 
modular polynomials. Furthermore, our algorithm also yields as a corollary a fast way to compute 
a random ^-isogeny of an elliptic curve over a finite field. 

The idea of our algorithm is as follows. Mestre's algorithm, Methode des graphes |Mes86j . uses the 
^th modular polynomial modulo p to navigate around the connected graph of supersingular elliptic 
curves over Fp2 in order to compute the number of edges (isogenics of degree t) between each node. 
From the graph, Mestre then obtains the Brandt matrix giving the action of the ^^'^ Hecke 
operator on modular forms of weight 2. In our algorithm we do the opposite: we compute the 
modular polynomial modulo p by computing all the isogenics of degree £ between supersingular 
curves modulo p via Vein's formulae. Specifically, for a given j-invariant, j (say), of a supersingular 
elliptic curve over Fp2, Algorithm 1 computes (j)g{x,j) modulo p by computing the I + \ distinct 
subgroups of order i and computing the j-invariants of the l + \ corresponding £-isogenous elliptic 
curves. Algorithm 2 then uses the connectedness of the graph of supersingular elliptic curves over 
Fp2 to move around the graph, calling Algorithm 1 for different values of j until enough information 
is obtained to compute (p£{x,y) modulo p via interpolation. 
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There are several interesting aspects to Algorithms 1 and 2. Algorithm 1 does not use the fac- 
torization of the ^-division polynomials to produce the subgroups of order i. Instead we generate 
independent ^-torsion points by picking random points with coordinates in a suitable extension of 
Fp and taking a scalar multiple which is the group order divided by i. This turns out to be more 
efficient than factoring the division polynomial for large 

Algorithm 2 computes (pe{x,y) modulo p by doing only computations with supersingular elliptic 
curves in characteristic p even though (f)i{x,y) is a general object giving information about isoge- 
nics between elliptic curves in characteristic and ordinary elliptic curves in characteristic p. The 
advantage that we gain by using supersingular elliptic curves is that we can show that the full 
^-torsion is defined over an extension of degree 0{i) of the base field Fp2, whereas in general the 
field of definition can be of degree as high as — 1. 

In this article we provide a running time analysis assuming fast multiplication implementation of 
field operations. But for small values of i fast multiplication is usually not used in practice, thus 
wc also give the running time (without the analysis) assuming a naive implementation of field 
operations. 



2. Local computation of 

The key ingredient of the algorithm is the computation of the univariate polynomial (pi(x,j) mod- 
ulo a prime p given a j-invariant j. We describe the method to do this here. 

Algorithm 1 

Input: Two distinct primes p and i, and j the j-invariant of a supersingular elliptic curve E over 
a finite field ¥q of degree at most 2 over a prime field of characteristic p. 
Output: The polynomial ^e{x, j) = H^, ^.i^ogenous to " ^'(^0) ^ Fp2 [x] . 

Step 1 Find the generators P and Q of E[i]: 

(a) Let n be such that ¥q{E[e]) C F^n. 

(b) Let S = '^E{¥qn ) , the number of F^n rational points on E. 

(c) Set s = S/i'^, where is the largest power of i that divides S (note A; > 2). 

(d) Pick two points P and Q at random from E[£]: 

(i) Pick two points U, V at random from E{¥qn). 

(ii) Set P' = sU and Q' = sV, if either P' or Q' equals O then repeat step (i). 

(iii) Find the smallest ii,i2 such that f ip' 7^ O and I'^Q' 7^ O but i^^+^P' = O and 

(iv) Set P = t^P' and Q = t^Q'. 

(e) Using Shanks's Baby-stcps-Giant-stcps algorithm check if Q belongs to the group gen- 
erated by P- If so repeat step (d). 

Step 2 Find the j-invariants j'l, • • • ,je+i in Fp2 of the £ + 1 elliptic curves that are ^-isogenous to 
E. 

(a) Let Gi = {Q) and d+i = {P + {i - l)Q) for 1 < i < £. 

(b) For each i, 1 < i < i + 1 compute the j-invariant of the elliptic curve E/Gi using Vein's 
formulas. 

Step 3 Output (t)e{x,j) = UiKiKe+ii^ - Ji)- 
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Remark 2.1. In step (le), one could alternately use the Weil pairing to check whether P and Q 
generate the ^-torsion. Doing so, however, does not lead to an asymptotic improvement in the 
running time of the algorithm. 

The following lemma gives the possibilities for the value of n in Step (la). We prove the following 
result for all elliptic curves not just supersingular ones. 

Lemma 2.2. Let E/¥q be an elliptic curve, and let i be a prime not equal to the characteristic of 
¥q. Then E[£] C E(¥gn) where n is a divisor of either £(£ — 1) or f?' — 1. 

Proof : The Weil-pairing tells us that if E[e\ C F^^n then p.^ C Fgn f p?il86] Corollary 8.1.1). We 
seek, however, an upper bound on n, to do this we use the Galois representation coming from the 
^-division points of E. Indeed, we have an injective group homomorphism f |Sil86j Chapter III, §7) 

pe : Gali¥g{E[£])/¥g) ^ Aut{E[i]) ^ GL2(F,). 

The Galois group Gal{¥ q{E[i])/¥q) is cyclic, thus by pi the possibilities for Gal{¥ q{E[£])/¥q) are 
limited to cyclic subgroups of GL2(F^). In other words, we are interested in the orders of the 
elements in GL2(F^). The elements of Gh2{¥£) are conjugate to one of the following types of 
matrices: 

(o °)'(o «)>fora,/3EF^a//3, 

or those corresponding to multiplication by an element of F^a on the 2-dimensional F^ vector space 
¥p. It is easy to see that the orders of these elements all divide £{£ — 1) or £^ — 1. Thus the de- 
gree of the field extension containing the ^-torsion points on E must divide either £{£—!) or £^ — 1. □ 

We will try step (1) with n = — 1, if steps (Id - le) do not succeed for some K (a constant) 
many trials, we repeat with n = £{£ — 1). The analysis that follows shows that a sufficiently large 
constant K will work. 

For step (lb) we do not need a point counting algorithm to determine S. Since £' is a supersingular 
elliptic curve, we have the following choices for the trace of Frobenius Oq: 

_ Jo if£;isoverFp 
[0,±p, ±2p if E is over Fp2 . 

Not all the possibilities can occur for certain primes, but we will not use this fact here (see |Sch87j ). 
If the curve is over Fp2 we can determine probabilistically the value of Og as follows. Pick a point 
P at random from E{¥q) and check if (g -|- 1 -|- aq)P = O. Since the pairwise gcd's of the possible 
group orders divide 4, with high probability only the correct value of Oq will annihilate the point. 
Thus in 0(log2+°(^) q) time we can determine with high probability the correct value of aq. Once 
we know the correct trace a^, we can find the roots (in Q), vr and vf, of the characteristic polynomial 
of the Frobenius cfP' — aq(f) + q. Then the number of points lying on E over the field ¥qn is given by 
-I- 1 — vr" — Tf*^, this gives us S. 

Note: We could have used a deterministic point counting algorithm to find ji£'(Fg) but this would 
have cost Oilo^ q) field operations. 

A lower bound on the probability that step (Id) succeeds is given by the following lemma whose 
proof is straightforward. 
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Lemma 2.3. For a random choice of the points U and V in step (Id i) the probability that step 
(Id ii) succeeds is at least 




At the end of step (Id) we have two random ^-torsion points of E namely, P and Q. The probabihty 
that Q belongs to the cyclic group generated by P is ^ = |. Thus with high probability we will 
find in step (le) two generators for E[i]. 

Lemma 2.4. The expected running time of Step 1 is 0{£^^°^^^ log^^"*-^^ q) and 0{i^ log^ q) if fast 
multiplication is not used. 

Proof : The finite field F^n can be constructed by picking an irreducible polynomial of degree 
n. A randomized method that requires on average 0((n^ log n + nlogg) log nloglogn) operations 
over ¥g is given in |Sho94j . Thus the field can be constructed in 0(£''+°(^) log^+"(^) q) time since 
n < l"^. Step (Id) requires picking a random point on E. We can do this by picking a random 
element in F^n treating it as the x-coordinate of a point on E and solving the resulting quadratic 
equation for the y-coordinate. Choosing a random element in ¥gn can be done in 0{£'^ logg) time. 
Solving the quadratic equation can be done probabilistically in 0{£'^ logq) field operations. Thus 
to pick a point on E can be done in 0(1^'^°^'^'^ log^"*""^^^ q) time. The computation in steps (Id i - 
iv) computes integer multiples of a point on the elliptic curve, where the integer is at most g", and 
this can be done in 0(£4+°(i) log2+°(i) q) time using the repeated squaring method and fast multi- 
plication. Shanks's Baby-steps-giant-steps algorithm for a cyclic group G requires 0(y^|G|) group 
operations. Thus step (le) runs in time 0(^2+°(^) log^'*'"*-"'^^ q), since the group is cyclic of order i. □ 

Let C be a subgroup of E, Velu in |Vel71j gives explicit formulas for determining the equation of 
the isogeny E — > E/C and the Weierstrass equation of the curve E/C . We give the formulas when 
£ is an odd prime. Let E is given by the equation 

+ aixy + osy = + 02^^ + + a^. 

We define the following two functions in ¥q{E) for Q = {x,y) a point on E — {O} define 

g^{Q) = 3x^ + 2a2X + 04 — aiy 
g^{Q) = -2y - aix - as, 

and set 

t{Q) = 2g^Q)-a,gy{Q) 
u{Q) = {gy{Q)f 

t = *(^) 

Qe(c-{ci}) 

w= Yl HQ) + x{Q)tm- 

Qe(c-{o}) 

Then the curve E/C is given by the equation 

+ AiXY + A3Y = X^ + A2X^ + AaX + Aq 
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where 



Ai = ai,A2 = 02, ^3 = as, 

A4 = a4 — 5t, Aq = uq — {af + Aa2)t — 7w. 

From the Weierstrass equation of E/C we can easily determine the j-invariant of E/C . It is clear 
that this procedure can be done using 0{t} elliptic curve operations for each of the groups G,, 
1 < i < (. + 1. Thus step 2 can be done in 0(£^+°(^) log^"*""^^^ q) time steps. Step 3 requires only 
0{tj field operations and so the running time of the algorithm is dominated by the running time 
of steps 1 and 2. Note that the polynomial obtained at the end of Step 3 (t)^{x^j) has coefficients 
in Fp2[x] since all the curves ^-isogenous to E are supersingular and hence their j-invariants belong 
to Fp2. In summary, we have the following: 

Theorem 2.5. Algorithm 1 computes (/)i{x,j) G Fp2[x], in fact, the list of roots of (j)g{x,j), and 
has an expected running time of 0{i^'^°^^^ log^^°^^^ q) and 0{£^\og^ q) without fast multiplication. 

For our application of Algorithm 1 we will need the dependence of the running time in terms of 
the quantity n. We make the dependence explicit in the next theorem. 

Theorem 2.6. With notation as above, Algorithm 1 computes (j)£{x,j) G Fp2[3;] together with the 
list of its roots and has an expected running time of 0{n'^~^°^^^ log^^"*-^-* q + v^n^^"^^-* log"*^"*""^"^^ q + 
£2^i+o{i) logfji). If fast multiplication is not used then Algorithm 1 has an expected running time of 
0{n^ log^ q + \/In^ log^ q + fn'^ log^ q). 

In the case of ordinary elliptic curve, step (1) of Algorithm 1 can still be used, once the number of 
points on E/¥q has been determined, by Lemma 12.21 the degree of the extension, n, is still Oi/^). 
This leads to the following two results: 

Corollary 2.7. // E/¥q is an elliptic curve, we can pick a random l-torsion point on E(¥q) in 
time 0(£^+°(i) log^+°(^) q + log^+°(^) q) and 0{e^ log^ q + log* q) without fast multiplication. 

Corollary 2.8. If E/¥q is an elliptic curve, we can construct a random l-isogenous curve to E in 
time 0(£^+°(^) log^+°(^) q + log'^+°(^^ q) and 0{i^ log^ q + log* q) without fast multiplication. 

Factoring a degree d polynomial over a finite field F^n requires 0{d^{nlogq)) operations over F^n. 
To factor the ^-division polynomial we need an extension of degree roughly i'^. Thus, if we were to 
factor the ^-division polynomial to generate the isogeny, we would need 0{£^ logq) operations over 
a field of degree i"^ over F^ which translates to 0{i^~^"^^^ log^"*""*-^-* q) bit operations even with fast 
multiplication. 

3. Computing y) mod ^» 
In characteristic p > 2 there are exactly 

S{p) 

supersingular j-invariants where 

ep = 0,1,1,2 if p = 1,5,7, 11 mod 12. 

In this section we provide an algorithm for computing (pe{x, y) mod p provided S{p) > £ + 1. The 
description of the algorithm follows: 

Algorithm 2 

Input: Two distinct primes £ and p with S{p) > i + 1. 
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Output: The polynomial (l)i{x,y) £ ¥p[x,y]. 



(1) Find the smallest (in absolute value) discriminant D < such that (-p) = —1. 

(2) Compute the Hilbert Class polynomial Hd{x) mod p. 

(3) Let jo be a root of Hd{x) in Fp2. 

(4) Set i = 0. 

(5) Compute (pi = (j)i{x,ji) G ¥p2 using Algorithm 1. 

(6) Let ji+i be a root of (pk for < i which is not one of jo, ■ ■ ■ ,ji. 

(7) If i < i then set i = i + 1 and repeat Step 5. 

(8) Writing (j)(^{x,y) = x^^^ + YliQ<k<iPk{y)x^ ■, we have systems of equations of the form 
Pk{ji) = Vki for < A;, i < i. Solve these equations for each Pk{y), < k < £. 

(9) Output (piix^y) e Vp[x,y]. 

We argue that the above algorithm is correct and analyze the running time. For step 1, we note that 

if p = 3 mod 4, then D = —4 works. Otherwise, —1 is a quadratic residue and writing (without 

loss of generality) D as —4d, we are looking for the smallest d such that (p) = ~1- A theorem 

1 

of Burgess ( |Bur62j ) tells us that d <^ p*^, and under the assumption of GRH the estimate 
of Ankeny ( |Ank52j ) gives d <C log^p. Computing Hd{x) mod p can be done in 0((i^(log d)^) 
time (LL9n] §5.10. Thus step 2 requires 0{y/plog^ p) time, and under the assumption of GRH 
requires 0(log^p(loglogp)^) time. Since (-p) = —1 all the roots of Hd{x) are supersingular j- 
invariants in characteristic p. H£){x) is a polynomial of degree h{y/—D), the class number of the 
order of discriminant D, and this is <C |Z)|2'^'^. Finding a root of Ho{x) £ ¥p2 can be done in 
0{d^+^ log2+"(^) p) time using probabilistic factoring algorithms, where d = \D\. The graph with 
supersingular j-invariants over charactertistic p as vertices and ^-isogenics as edges is connected (see 
|Mes86j ). consequently, we will always find a j-invariant in step 6 that is not one of jo, • • • ,ji. Thus 
the loop in steps (5) • • • (7) is executed exactly £ + 1 times under the assumption that S{p) > £ + 1. 
Even though Algorithm 1 requires O(^^log^g) time^ in the worst case, we will argue that, in fact, 
for all of the iterations of the loop it actually runs in 0{£^ log^ q) time. 

Lemma 3.1. Let E be a supersingular elliptic curve defined over¥p2. Then the extension degree 

[¥p2{E[£]) : ¥p2] 

divides 6{£ — 1). 

Proof : Let E/¥p2 be a supersingular curve and let t be the trace of Frobenius. Then the Frobenius 
map (j) satisfies 

- t0 + p2 = 

with t = or zizp or ±2^. Suppose t = ±2p, then the Frobenius acts as multiplication by itp on the 
curve E. Thus (p^~^ acts trivially on E[£], and the ^-torsion points are defined over an extension 
of degree dividing £ — 1. If t = 0, then (p'^ = —p^ and so acts trivially on the ^-torsion. 

Thus E[£] is defined over an extension of degree dividing 2{£ — 1). If t = itp, then cp^ = itp^ and 
consequently acts trivially on the ^-torsion of the curve. Thus the f-torsion is defined over 

an extension of degree dividing 2>{£ — 1). Thus in all cases the ^-torsion of the curve is defined over 
an extension of degree dividing 6(£ — 1). □ 



We use the soft-Oh O notation when we ignore factors of the form log^ or log log p. 
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Thus, the loop in steps (5) • • • (7), Algorithm 1 can be run with the quantity n = 6{£ — 1). For this 
value of n Algorithm 1 runs in expected time 0{£^~^°^^^ log^"*""^^^ p) and so the loop runs in expected 
time 0(£4+"(i) log2+°(^)p). 

Writing the modular polynomial (f)g{x,y) as x^~^^ + 'Yl,Q<k<(Pk{y)^'^ we know that po(y) is monic 
of degree i+1 and deg(pfc(y)) < I for 1 < k < i. Thus at the end of the loop in steps (5) • • • (7) we 
have enough information to solve for the Pk{y) in step (8). We are solving £+1 systems of equations, 
each requiring an inversion of a matrix of size (^ + 1) x (£+ 1). This can be done in 0(^^ log^"*""^^^ p) 
time. Since the polynomial (j)e{x,y) mod p is the reduction of the classical modular polynomial, a 
polynomial with integer coefficients, the polynomial that we compute has coefficients in ¥p. Thus 
we have proved the following theorem: 

Theorem 3.2. Given £ and p distinct primes such that S{p) > £ + I, Algorithm 2 computes 
(f)£{x,y) G Fp[x,y] in expected time ©(^^"'""^''^^ log^"*""^^^ p + log^ploglogp) under the assumption of 
GRH and in 0{£^ log^p) time without fast multiplication. 

Hence, we can compute (j)^{x,y) modulo a prime p in 0(^^log^p + log^p) time ii p > 12£ + 13. If 
p < 12£ + 13, we could still use the algorithm with ordinary elliptic curves and this would lead to a 
running time with the dependence on £ being £^. Furthermore, we would not need the GRH since 
it was needed only to determine a supersingular curve in characteristic p. However, this is not very 
efficient. 

Rem,ark 3.3. If one is allowed to pick the prime p, such as would be the case if we arc computing (pn 
over the integers using the Chinese Remainder Theorem combined with this method, then one can 
eliminate the assumption of GRH in the above theorem. For example, for primes p = 3 mod 4 the 
j-invariant 1728 is supersingular. Thus in step (3) of Algorithm 2, we can start with jo = 1728 for 
any such prime. Hence we do not need the GRH to bound D in the analysis of the running time 
of the algorithm. 
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Appendix 

Elkies in |Elk98j (see §3 of that paper) claims that the usual method of computing the i-th modular 
polynomial runs in time 0(£^"^'^). However, there is a subtle error in the analysis. We argue that 
in fact, the running time of this algorithm is 0(^2+'^). The first stage of the algorithm invovles 
computing the first l"^ + 0{£) Fourier coefficients of the powers of the j-function, namely j, • • • , j^- 
This (as Elkies points out) can be done in 0(^'^+'^) arithmetic operations. The problem comes when 
we study the running time in terms of bit- operations. To analyze this we need to study the bit-sizes 
of the numbers that are handled by the algorithm. While it is true that the n-th Fourier coefficient 
of j grows as e^^^^ , we need to also compute the Fourier coefficients of powers of j and they grow 
faster (essentially because they have a higher order pole at oo). In |Mah74j an upper bound of the 
form exp(47r(y^ (n -|- k)k) is proven for the n-th Fourier coefficient of j'^. We show that the upper 
bound is quite close to the true magnitude of the Fourier coefficients below. Let c(n) denote the 
n-th Fourier coefficient of the j-function. It is well known that c(n) are all positive integers and 
that (see |Pet32p 

g47rv^ 

The n-th Fourier coefficient of j'^ is given by 

c(ai)c(a2) • • • c(afc). 

ai-l-ct2-l l-afe=n 

Clearly, there is at least one partition of n (into k parts) where each of the parts, Oj, are > 
where c > 1 is a constant. The asymptotic formula for c(n) gives us that 

J2 c(ai)c(a2)---c(afc) > (^eV^~Oi^°sn)y ^ ^n{V^) ^ 

21+12-1 |-afc=n 

as long as k and n vary such that the ratio n/k goes to infinity. Thus a lower bound for the rate 
of growth of the n-th Fourier coefficient of j'^(z) is e^^^^^\ Hence to compute the first i'^ + 0{£) 
Fourier coefficients of the powers j'^, for 1 < k < i, the bit length of the numbers involved is 0{iVi) 
rather than the estimate O(^log^) used in |Elk98j . Thus this stage of the algorithm already requires 
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0(^2"'"'^) time using fast multiplication. The dependence of the running time on £ for our method, 
on the other hand, is 0{£^^'') if fast multiplication is used. The error in |Elk98j stems from the 
fact that a bound on the size of the coefficients of (j)i{x,y) was used to bound the bit sizes. These 
coefficients are somewhat smaller, their absolute value being bounded by e^(^'°g^) (see |(^ohS4) ). 
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